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Continuous Markovian model for Levy random walks 
with superdiffusive and superballistic regimes 
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Abstract. We consider a previously devised model describing Levy random walks (Phys. Rev E 79, 011110; 
80, 031148, (2009)). It is demonstrated numerically that the given model describes Levy random walks with 
superdiffusive, ballistic, as well as superballistic dynamics. Previously only the superdiffusive regime has 
been analyzed. In this model the walker velocity is governed by a nonlinear Langevin equation. Analyzing 
the crossover from small to large time scales we find the time scales on which the velocity correlations 
decay and the walker motion essentially exhibits Levy statistics. Our analysis is based on the analysis of 
the geometric means of walker displacements and allows us to tackle probability density functions with 
power-law tails and, correspondingly, divergent moments. 



PACS. 05.40.Fb Random walks and Levy nights 
processes - 05.10.Gg Stochastic analysis methods 



1 Introduction 

The Levy-type behavior of anomalous Brownian motion 
has been found in a large variety of systems in nature 
(see, e.g., [U^EJEUS])- It is characterized by the scaling 
law \5t\ oc t& of the particle displacement \Sr\ during the 
time interval t with an exponent j3 > 1/2. The cases of 
1/2 < f) < 1 and (3 > 1 are usually referred to as the 
superdiffusive and superballistic regimes of the Levy-type 
dynamics, respectively. 

The continuous time random walk model, including its 
decoupled and coupled implementations, is now a widely 
used approach to describe Levy-type stochastic processes 
(see, [SJE] and references therein) . Within the classical for- 
mulation Levy flights are modeled as a sequence of random 
independent steps whose lengths are distributed accord- 
ing to a power-law probability density with diverging mo- 
ments. Unfortunately, this model does not generate con- 
tinuous trajectories. Indeed, by one step a Levy walker 
can make a very long jump and such events in spite of 
being rather rare contribute substantially to the walker 
displacements. This feature is responsible for considerable 
problems met in describing Levy-type processes in hetero- 
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geneous media and systems with boundaries. The coupled 
continuous time random walks [B] , also referred to as Levy 
random walks, do admit a representation of the walker mo- 
tion in terms of continuous trajectories; this model just 
mimics the stochastic process at hand as a sequence of 
random independent steps whose length and duration are, 
in turn, random variables. The latter enables one to in- 
troduce a finite velocity of the walker motion during one 
step by connecting its initial and terminal points with a 
straight line and considering the walker motion along it 
uniform. Unfortunately, the resulting stochastic process 
looses the Markov property in the strict sense. 

There have been several attempts to develop a descrip- 
tion of Levy-type random walks starting from a certain 
Markovian continuous process with, may be, nonlinear 
(multiplicative) Langevin forces. In particular, it is a de- 
scription of the decoupled continuous time random walks 
based on the corresponding Langevin equation [7], the 
combination of a conventional stochastic process governed 
by a linear Langevin equation with additive noise where, 
however, time itself is a random nondecreasing variable 
8 , models bases on stochastic differential equations with 
non-Gaussian Levy white noise [9,10 , and the semiclas- 
sical Fokker-Planck equation describing quantum interac- 
tion of atoms with laser beams pTl[T2] , as well as models 
using the Boltzmann equation with scattering character- 
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ized by power- law distribution [T31[T4"1[T5] . There is also 
an approach to describe anomalous diffusion based on the 
nonlinear Fokkcr-Planck equation derived using the gen- 
eralized entropy |16| . 

In systems far from thermodynamics equilibrium and 
complex systems nonlinear stochastic processes are widely 
met and have been studied in connection with the appear- 
ance of power-law distributions and long-time tails of time 
correlation functions (see, e.g., [TTlfTB] as well as [T91I50"] 
and references therein). 

The generalized Cauchy stochastic process, which is 
governed by the following Langevin equation with additive 
and multiplicative noise [20] 

T-£ = -[\ + ^U(t)]*V + y/%fVaZ a (t) (la) 

at 

has been employed to study stochastic behavior of various 
nonequilibrium systems, in particular, lasers [T7], on-off 
intermittency [21], economic activity [22], passive scalar 
field advected by fluid [53], etc. Here v is a certain ran- 
dom variable, r is a characteristic "microscopic" time scale 
of its dynamics, A is a friction constant, ^ m (t) and t; a (t) 
are independent components of white noise of unit ampli- 
tude, the parameter v a quantifies the intensity of additive 
noise, and the asterisk indicates that equation dial) is writ- 
ten in the Hanggi-Klimontovich form. The coefficients \f~2 
are introduced into (fla|) for the sake of convenience. Equa- 
tion (|la.p is equivalent to the following Langevin equation 
with a single white noise source £(£) and the amplitude 
depending nonlinearly on the random variable v |20] 



dv 
T—- = —Av 
dt 



2t(i 



2 W 



(lb) 



because both of them are described by the same Fokker- 
Planck equation for the probability density P(i>, t) 



dP 
"dt 



d 


8P 


dv 


{< + v 2 )-^ + \vP 



(2) 



Appealing to qualitative arguments and a numerical ex- 
ample Sakaguchi [53] demonstrated that the generalized 
Cauchy process (fTa|) can give rise to Levy flights in the 
space x on time scales t 3> t if the random variable v is 
regarded as the velocity v = dx/dt of a Brownian particle. 
Previously 25,26 we analyzed in detail such random 
walks governed by the Langevin equation (|lb[) when the 
induced Levy-type process is superdiffusive. In particular, 
it has been demonstrated [55] that on time scales t S> r, 
related with the crossover to Levy-type behavior, the scal- 
ing law 5x ex t 13 with the exponent 1/2 < /3 < 1 is due 
to the impact of extremely large fluctuations in the ve- 
locity pattern {v(t)}. Roughly speaking, the peak velocity 
{fmax(i)} times the duration r, where v m&x is the maximal 
value of the velocity attained during the time interval t, is 
the main contribution to the particle displacement 8x(t). 
As a consequence, the estimate Sx ~ v maK T holds and 
the characteristic scaling law is related to the statistics of 
the occurrence of extreme fluctuations in the observation 
interval t. 



It should be noted that the constructed random pro- 
cess neither belongs to the decoupled nor to coupled con- 
tinuous time random walks. Indeed, on one hand side, this 
Brownian particle moves continuously and, on the other 
hand, its motion with the "elementary" time fragment t 
is not uniform but extremely heterogeneous. The latter 
is actually the mechanism preserving the Markov prop- 
erty of such random walks. Naturally on time scales about 
t and the corresponding spatial scales £ — v a r random 
walks governed by the Langevin equation (|lb[) are ballis- 
tic, Sx(t) oc t, because of the strong velocity correlations 
on such scales. In our paper [5(5] we have developed a spe- 
cial singular perturbation technique and proved, for the 
superdiffusive regime, that the stochastic process ATE]) in- 
deed describes the Levy-type random walks on time scales 
t 3> r. At the critical value A c separating the superdiffu- 
sive and superballistic regimes the constructed perturba- 
tion technique exhibits some singularities at intermediate 
steps, but the final results turn out to be free from any 
singularity at A = A c . This problem seems to be related 
to the well-known conflict with dynamics in constructing 
the superballistic Levy flights, extensively discussed in the 
review [27] , 



The purpose of the present paper is, first, to demon- 
strate that the developed model 25 , 26 not only describes 
the superdiffusive Levy-type random walks but also su- 
perballistic ones. A second aim is to analyze the crossover 
between the time scales t ~ t and t ^$> t in the behavior of 
the walker displacements. The latter result will clarify the 
behaviour at the time scales for which velocity correlations 
become insignificant and Levy-type dynamics arises. 



It should be noted that superballistic transport is ex- 
pected to arise in such systems as atmospheric turbu- 
lence [551155] , two-dimensional turbulence [3U] , a collision- 
less plasma driven in magnetic field reversal with an elec- 
tric field in the presence of magnetic turbulence |31j . It 
also has been observed experimentally in a fully devel- 
oped turbulent laboratory water flow [35] as well as in 
a one-dimensional expansion of cold atomic clouds in a 
magneto-optical trap under the impact of laser beams |33) . 
In spite of these objects, strictly speaking, do not belong 
to the systems under consideration, it was found that the 
notion of Levy-type random walks [291130] or stochastic 
differential equations [33J gives rise to the adequate de- 
scription of these superballistic phenomena. Large-scale 
animal movements in searching resource patches, mates, 
etc. or dispersal also exhibit the Levy- type pattern (for a 
review see 34,35,36 ) and there are conditions when the 
optimal random search becomes of the ballistic type as a 
limit case of Levy- type motion [37,38 . Levy type behavior 
is relevant to stock indices, too, and the indices matching 
the superballistic Levy flights were found, in particular, 
for the Brazil BVSP index and the Argentine MERV in- 
dex [39], for foreign exchange rate in Brazil and South 
Korea [4TJU4"T] , and the absolute value of the returns in the 
Indian market indices [35] 
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2 Model 

Let us restrict our consideration to a dimensionless de- 
scription of the stochastic processes at hand and, thereby, 
measure time in units of r, i.e. t —¥ rt and spatial vari- 
ables in units of £, i.e. x% — > ixi. In these units the con- 
tinuous Markovian random walks under consideration are 
governed by the following equations 



dxi 

~~dt 
dvj 
~dt 



(n + a)v t + v/2(l + v 2 ) * &(£) . 



(3) 

(4) 



Here r = {x\, . . . ,x n } is the current position of a ran- 
dom walker in the Euclidean n-dimensional space R" and 
v = {vi , . . . , v n } is the velocity of its motion in this phase 
space, {£i(£), • • • £n(*)} is the collection of mutually inde- 
pendent components of Gaussian white noise, 

(&(t)}=0, and (£i(t)&(t'))=6 ii >6(t-t'), (5) 

and, as before, the asterisk indicates that equation (j4]) is 
written in the Hanggi-Klimontovich form. The parameter 
a is assumed to be taken from the interval 



< a < 2. 



(6) 



Previously [55] dealing with this model for 1 < a < 2 
we demonstrated that in the limit t ^> 1 the characteris- 
tic amplitude |<5r| of the walker displacements scales with 
time t as | <5r(i) | ex t 1 /", matching the superdiffusive Levy 
random walks. Besides, the probability density P(5r, t) of 
the walker displacement Sr £ R™ is found to possess the re- 
quired asymptotics P(Sr,t) ex |Jr|~' n+Q! ) when \5r\ —> oo. 
These expressions stem directly from the obtained depen- 
dence of the generating function for t 3> 1 



G(k,t) = <e lkr ) = exp{-A mJn t|k| a } ) 



(7) 



where yl m i n > is the minimal eigenvalue of the cor- 
responding Fokkcr-Planck equation. In the ID-case it is 
given by the expression 



r 



A 



2-a 



mm — /q,\ 

2^r[-)r(a) 



(8) 



where r(. . . ) is the gamma function. 

As noted before all these final results do not exhibit 
singularities when a — ► 1 whereas the constructed pertur- 
bation technique does at its intermediate steps. This fea- 
ture poses a question as to whether the stated model ([3]), 
Q does describe the Levy random walks in a wider region 
of parameters including the whole interval (JB]) and com- 
prising superdiffusive as well as superballistic dynamics. 

In the next section this hypothesis will be justified nu- 
merically; here, we present some heuristic arguments for 
the absence of any singularity in the behavior of the given 
system when the parameter a is taken from the interval 
< a < 1. 



Canonical form of equation Q 

Dealing with the one dimensional stochastic processes gov- 
erned by the Langevin equation (01) for n — 1 let us make 
use of the nonlinear transformation 1 20 



V = sinh(w) 



(9) 



reducing equation (jj) to the following Langevin equation 
with additive noise 



da 
~~dt 



-^ = -atanh(w) + \/2f(t) 



(10) 



In obtaining expression (|10j) we, first, converted equa- 
tion (J4]) into the Stratonovich form and only then per- 
formed the transformation ©. According to Konno and 
Watanabe [3D] expression (flU)) can be regarded as the 
canonical form of such stochastic differential equations 
with multiplicative noise. 

From the Fokker-Planck equation 



dP 
~dt 



d_ 

Oil 



dP 
du 



a tanh(u)P 



corresponding to the Langevin equation (|10p we obtain 
the steady state distribution of the random variable u 



Pc(u) 



A r 



[cosh(u)] a ' 



(11) 



where A c is the normalization constant. In particular, for 
a ~ 1 the constant A c ~ 1 too and for a <§C 1 the estimate 
A c = a/2 holds. From the Langevin equation (fTU|) we can 
estimate the time scale T c characterizing the formation of 
the steady state distribution pT|) as T c ~ 1 for a ~ 1 and 
T c - 1/a 2 for a < 1. 

We conclude that the statistical properties of the ran- 
dom variable u cannot exhibit any singularity for the pa- 
rameter a lying inside the interval under consideration, 
i.e. a € (0,2). The parameter a just affects the weight 
of extreme fluctuations with respect to their contribu- 
tion to the walker displacement, changing the exponent of 
its time dependence. Moreover, in studying these extreme 
fluctuations we can approximate the function tanh(u) by 
the signum function sign(x), i.e., make use of the replace- 
ment tanh(u) — > sign(x) and perform the transformations 
u — > u/a and t — > t/a 2 . This reduces equation (flU|) to the 
following equation 



du 

~di 



= -sign(u) + V2£(t) . 



(12) 



In deriving this equation the white noise transformation 
£(t) — > a£(t) caused by the given time renormalization and 
stemming from the induced transformation of the correla- 
tion function 



(mm) = s(t-t>) 



¥t/c 



a 2 5(t-t') 



has been taken into account. Equation (jT2|) does not con- 
tain any parameter, which seems to be the reason for the 
universal behavior of such systems. 
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Therefore, it is quite natural to expect that for < 
a < 1 model ([3]) , ((U) also describes the Levy-type random 
walks with, maybe, the superballistic dynamics provided 
the expression 5x{t) ex t x l a holds. The singularities arising 
at intermediate steps of the perturbation technique seem 
to be caused by anomalous behavior of integrals similar 
to 



and for the 2D-case 



P(x, t) sin(fca:) dx oc 



sin(fcx) 



ll + Q 



dx 



as k — > 0. For a < 1 this integral also leads to nonanalytic 
dependencies on the wave number k and must vanish for 
the symmetry reasons only. 



3 Numerical results 

In order to justify the fact that the model (J3J), (J3J) does 
describe the Levy-type random walks on time scales t > 1 
for < a < 2 the system of equations ©, (01 has been 
solved numerically. Then the obtained results were used 
to calculate the geometric mean value 



ln[F g (t)]=(ln(|<5r|)> 



(13) 



of the walker's displacement <5r during the time interval 
t as well as the probability density function P(\Sr\,t). It 
should be pointed out that logarithmic mean values simi- 
lar to (|13[) form a rather efficient tool for analyzing statisti- 
cal data with power-law distributions possessing diverging 
moments 20 . It is due to the geometric mean converg- 
ing for any power-law distribution that can be normalized 
to unity and, thus, enables us to analyze superdiffusive 
and superballistic regimes with any possible value of the 
exponent a € (0, 2) using the same quantities. The geo- 
metric mean value introduced via expression (| 13[) can be 
also treated as the limit case of constructing the Hurst 
exponent H(e) 



(\5r\ 



l/« K t H{e) 



for 



-0. 



In addition, when the random variables Sr(t) are charac- 
terized by an autonomous probability density V (|<5r|/£' 3 ) 
expression (fT3|) reads 



ln[r„(t)] = jHQViOdC + PHt), 
o 

which allows us to separate the time dependence of the 
scaling low |<5r(£)| oc t@ from the corresponding numerical 
coefficients. 

If the generating function G(k, i) for the quantities 
{Sr} t is given by expression Q formula (|L?)) yields 



1 {t) = A g .t l ' a , 



(14) 



as shown in Appendix O In particular, we obtain in the 
ID-case the coefficient 



A g \ n=1 = A]i" n ■ exp <J 7 



(I -a) 



(15) 



A n 



2^-exp 



7- 



(1 



(16) 



Here 7 w 0.5772 is the Euler-Mascheroni constant. 

If the random variable Sr obeys the Levy statistics, 
then the asymptotics of the probability density P(\8r\,t) 
is also determined by equation (JT]) for the generating func- 
tion G(k, t). Moreover, provided the walker displacement 
Sr is measured in units of the corresponding geometric 
mean r g (t), i.e. in terms of p = \5r\/r g (t), its probability 
density is characterized by the asymptotics (Appendix [Bj) 



HP) 



nl + Q 



as p — > 00 , 



where in the ID-case the coefficient reads 
w a \ n=1 = ~ sin (^) r (1 + a) e^" 1 ) 



whereas in the 2D-case 



VV a \ n= 2 



2 . /na\ 

— sin — — 
n V 2 / 



a 

/U+ 2 



P 7(a-1) 



(17) 



(18) 



(19) 



Expressions (fT4)l and (fT7| are the key results for the nu- 
merical verification of the Levy-type dynamics of the ran- 
dom walks under consideration. Namely, analyzing the nu- 
merical simulation data it is possible to verify whether the 
constructed time dependence of the geometric mean r g (t) 
and the probability density P(r, t) actually do obey ex- 
pressions (HU) and (1T7|) . 

For the ID-case the system of equations (O, (@| has 
been solved numerically employing the stochastic Runge- 
Kutta method SRI2W1 approximating one-dimensional 
Wiener process Ito equations with strong order (4.0,1.5) 
43]. The time step dt = 0.05 for integration yields results, 
which were stable with respect to a decrease of dt. The 
integration time was T = 10 8 . For the 2D-case this sys- 
tem was solved using the stochastic Runge-Kutta method 
RI1 approximating multidimensional Wiener process Ito 
equations with weak order (3,2) [33]. The integration time 
step was chosen as dt = 0.02 yielding stable results. The 
integration time was also T = 10 8 . 

Figure Q] depicts the obtained results for the ID and 
2D models. The relevant analytic time dependence of the 
geometric mean x g (t), in the ID-case specified by expres- 
sions (0), ([T3|). and ([151), is also exhibited in Fig. [T] In 
particular, for a — 1.5, 1.0 and 0.5 we obtain A g w 1.46, 
1, and 0.130, respectively. As seen in Fig.[T]the correspond- 
ing analytic dependence fits the numerical data quite well 
if the time interval t exceeds a value of about t ~ 100. For 
time scale t ~ 1 at which the velocity correlations should 
be significant and, consequently, the walker behaviour has 
to be ballistic, the effective exponent of the curve x g (t) 
is about unity as expected. Thus, the constructed curves 
actually describe the crossover from the ballistic motion 
on time scales t ~ 1 to the superdiffusive or superballistic 
regimes of the Levy-type random walks for t > 100 for all 
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Fig. 1. The geometric mean of the walker displacements versus the observation time interval and the probability density of 
the walker displacements. In simulating ID and 2D models the integration time steps was set equal to dt = 0.05 and dt = 0.02, 
respectively, the integration time was T = 10 s . The other parameters are shown in the panels. 
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the analyzed values of the parameter a. In the 2D-case 
an analytic expression for the minimal eigenvalue yl m i n is 
not known. Therefore, we can compare only the slope of 
curves (in log-log scale) evaluated numerically and the cor- 
responding exponents of the theoretical dependence (IT4")) . 
As seen in Fig. Q] these values again fit each other when 
£ > 100 for the parameter a < 1 as well as a > 1. 

The numerical data for the probability density also 
exhibit the desired asymptotics for both the ID- and 2D- 
cases. According to expression (fT8|) in the ID-case the co- 
efficient w p takes the values zu p ~ 0.799, 0.300 for a = 1.5, 
0.5, respectively, and w p = 2/ir for a = 1. For the 2D-case, 
by virtue of ([T9"]l. it is 0.507, 0.277 and 1/2, respectively. 
The corresponding theoretical asymptotics of the proba- 
bility density V(p) are also shown in Fig. [TJ In particular, 
as seen in this figure the function V(p) tends to asymp- 
totics (fTTf from above for a = 1.5. 

For the parameter a — 1.5 the accuracy between the 
analytical and numerical results is expected and can be 
regarded as a verification of the used implementation of 
the employed numerical algorithms. At the same time, 
for a = 1.0 and 0.5 it is this accuracy that justifies the 
Levy-type behavior of the stochastic process governed by 
model ([3]), (HJ) for the parameter a lying in the interval 
< a < 1. 

It is necessary to note that the similarity of the plots 
for the ID and 2D-cases with respect to the correspond- 
ing exponents is a consequence of the fact that the friction 
coefficient (a + n) entering equation ((4]) depends on the 
space dimension n. So in the 2D-case the components £i(£) 
and £2 (t) of white noise individually undergo stronger dis- 
sipation and only their interplay leads to the same scaling 
law of the particle displacement magnitudes {|<5r| t }. The 
difference between the plots showing the time dependence 
of the geometric means x g (t) and r g (t) for a = 1 en- 
ables us to claim that the ballistic-like behavior of the 
geometric mean on time scales t < 1 and its ballistic 
regime for a = 1 must be characterized by the same ex- 
ponent but not the numerical cofactors. So their visible 
equality in the ID-case seems to be just coincidental. It 
also should be noted that in the 2D-case the probabil- 
ity density V{p) attains its maximum at a certain point 
Pmax > 0, whereas in the ID-case this maximum is lo- 
cated at the origin p = 0. In fact, even in the 2D-case 
the probability density P(Sr,t) — P(\5r\,t) of the random 
vector <5r gets its maximum at the origin too. However, 
the depicted probability density V{p) is proportional to 
the function P(\Sr\,t) with a coefficient depending on |<5r|, 
namely, V(p) ex |<5r|P(|<5r|, £), which causes the shift of its 
maximum from the origin. 

Finalizing the present section let us discuss ergodic 
properties of the investigated stochastic process. Ergodic- 
ity usually means that the ensemble average 




1(T 10" 10' 10' 10° 10° 10'" 
displacement Sx, in units of x g (t) 

Fig. 2. The two side probability density of the walker displace- 
ments. The data for the ID-model with the parameter a = 0.5 
matching the superballistic regime are presented. The integra- 
tion step was set equal to dt = 0.05 and the integration time 
was T = 10 8 . 

where P(A) is the probability density, and the time aver- 
age 



A:=- I A(t')dt' 



(A) := f AP(A)dA. 



of an observable A are equal in the infinite-time limit t —¥ 
00. The condition 

/(A-(A)) 2 \^0 as t^oo 

is widely used as a criterion for this equality 05]. A sys- 
tem obeying this criterion is said to be ergodic in the 
mean square sense. Therefore, the notion of ergodicity 
also includes the measure of convergence. This type of 
ergodicity can be broken when the parameter A is dis- 
tributed with a probability density possessing power-law 
tails which cause the divergence of its moments 12,46, 
147] . Moreover, studying the ergodicity one may focus not 
on a single trajectory of the system motion as a whole but 
treat it as a collection of various fragments. As a result, a 
system could be regarded as ergodic with respect to one 
variable and loose the ergodicity according to the proper- 
ties of another variable. Such a partial ergodicity is met 
in the system at hand. In fact, each of the plots in Fig. [1] 
has been built using numerical data generated during the 
system motion along a single sufficiently long trajectory. 
So the system ergodicity with respect to quantities such 
as the geometric mean of walker displacements and their 
distribution function is justified by the accuracy between 
the presented numerical data and the theoretical results 
related to the corresponding ensemble averages. As illus- 
trative example, Figure[2]depicts the two-sided probability 
density P(5x, t) of one-dimensional random walks for the 
parameter a = 0.5, i.e. for the superballistic regime. This 



Ihor Lubashevsky et al.: Continuous Markovian model for Levy random walks 



probability density, in contrast to the preceding construc- 
tions, takes into account also the direction of the walker 
displacements. As seen in this figure even in the case of the 
superballistic regime the probability density of the walker 
displacements constructed by averaging over time is highly 
symmetric. In contrast, under the same conditions such 
ergodicity does not hold with respect to any quantity de- 
pending on the sign of the walker position at the a;-axis. 
Indeed, the mean time t a spent by the random walker in- 
side a small neighborhood Q a of its initial point xq = 
whose size is about a can be estimated as 



t a 



X g (t) 



(It. 



(20) 



Estimate (|20p stems directly from the expression for the 
walker residence time t a inside the domain Q a 



t a = dt dx8[x — x{t)\ 



(21) 



for a given realization {x(t)} of the walker trajectory and 
the estimate of the probability density of finding the walker 
at the initial point x = 



g(o,t) = (S[x-x(t)])\ x=0 



i 



,(*) 



Integral (|2U)) converges at the upper boundary for the 
parameter a < 1 since in this case the geometric mean 
Xgit) oc t 1 /" grows faster than t. As a result, there is a 
certain probability that the walker never returns to its 
initial point and, thus, the walker will practically spend 
all the time in the region i>0ori<0ast- >oo. As 
the averages over ensemble are concerned, the properties 
of the given random walks should be symmetrical with 
respect to the reflection x i— > —x. 



4 Conclusion 

We have considered the stochastic process governed by 
equations ([3]), (jj) and numerically demonstrated that on 
time scales t ^$> 1 this process describes the Levy-type 
random walks with the superballistic, ballistic, and su- 
perdiffusive dynamics, depending on the value of the pa- 
rameter a, namely, for < a < 1, a = 1, and 1 < a < 2, 
respectively. In this way we have justified that the given 
model is a continuous Markovian implementation of the 
Levy flights or the Levy walks of all the main types that 
at the "microscopic" level admits the notion of continuous 
trajectories. 

To analyze the statistical properties of stochastic pro- 
cesses for which the moments of the probability density 
can diverge the notion of geometric mean was used and 
the walker displacements were measured in units of the 
corresponding geometric means. This enabled us to com- 
pare the numerical data and theoretical results not only 



for one-dimensional systems but also for two-dimensional 
models exemplifying the validity of the models for the case 
of multi-dimensional random walks. 

The analyzed crossover between the time scales t ~ 1 
and t > 1 has shown that the Levy-type behavior arises 
on scales t > 100. On scales t ~ 1, where the velocity cor- 
relations are significant, the random walks are ballistic. 
The direct simulation of the system motion along individ- 
ual sufficiently long trajectories has demonstrated the fact 
that the system at hand is partly ergodic, i.e. it is ergodic 
with respect to quantities such as the geometric mean of 
the walker displacements and their distribution, indepen- 
dently of whether the regime being superdiffusive or super- 
ballistic. In the latter case, however, there are quantities 
for which ergodicity is broken. At least the superballistic 
Levy-type random walks described by the given model are 
only partly ergodic. 

The authors are grateful to Hidetoshi Konno for his helpful 
comments and appreciate the financial support of the transre- 
gional research collaboration TRR 61 and the University of 
Munster as well as the partial support of the RFBR Grant 
09-01-00736. 



A Geometric mean and generating function 

Let us consider a random vector r of the n-dimensional 
space W l (n — 1,2,...) whose statistical properties are 
described by a given probability density P{v). Using this 
function we, first, construct the generating function G(k) 



G(k) := (e 4kr ) = I e 2kr P(r) dv , k e 



(22) 



which is the Fourier transform of the probability density 
and, then, the function 



G s (k) 



1 



^'n,s 



G(k') 



|k-k'| n - 



dk' . 



(23) 



where the exponent < s < 1, the coefficient 



2 S T 



^n,s = 7T 2 



r 



(24) 



and, as before, -T(. . .) is the gamma function. The substi- 
tution of (|22|) into (|23p yields the equality 



G s (k)= / e tkr P(r)dr 



1 



-■ -i i 



■ dq 



(25) 



A direct calculation of the last integral leads us to the 
expression 



G.,(k,= / ^P(r)dr. 



(26) 
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So, by definition, the mean value of the moment l/|r| s is 
1 \ f P(r) . 1 f G(k) 



Irl- a, 



\Un-a 



dk 



(27) 



by virtue of (J23J and (J2BJ. 

Let us introduce the symmetrized generating function 
g(k) by integrating the generating function G(k) over the 
n-dimensional sphere S nt k = {k : |k| = k} 



g(k) 



u n K 



L_^G(k)der, 

Sfc 



(28) 



where S n k n 1 is the surface area of the sphere §„,& and 
the coefficient 

n 

2tt2 



■5ji 



n? 



(29) 



In particular, for the random vector r with the spheri- 
cally symmetric distribution density P(r) and, thus, the 
spherically symmetric generating function G(k) the iden- 
tity g(k) = G(k) holds. In these terms expression (|2T| 
reads 



1 






k s - l g{k)dk 



*J n 



J /rj 



dk 



(30) 



Then using the series 



r 



|r|- s 


= 1- 


-sin 


\r\ + 0(s 2 ) 


fc s 


= l + slnk + 0(s 2 ) , 


2 s 


= l + ,sln2 + 0(s 2 ) , 


r (i) 


1 

s 


)-ll + 0(s*) 


77, — s\ 


-'(?) 


MGH 


2 j 



0(s 2 ) 



where, \P(. . .) is the digamma function and, as before, 7 ~ 
0.5772 is the Euler-Mascheroni constant, we expand the 
left and right sides of expression (I3TH) into series of s up 
to the first order in s. In this way for the geometric mean 
value r e of the random vector r defined as 



In [r g ] := (In |r[ 
we obtain the formula 



(31) 



ln[r g ] 



Hk)^dk- 

dk 



7 + 



V 



(?) 



$ 



(32) 
relating the geometric mean value r g and the generating 
function G(k) via g(k). In deriving expression (|32p the 
equality 

'/'' I -I =-7-21n2 



has been also taken into account. In particular, in the ID- 
case the last summand is absent and for the 2D-case where 
#-(n/2) = #(1) = -7 its value is [#(1) - #(l/2)]/2 = 
In 2. It should be noted that expression (|3"2j) can be also 
regarded as a specific implementation of the Hadamard 
type fractional calculus (see, e.g., |38]). 

For the generating function G(k) = exp (— oj|k| Q ), where 
u> > and a > are some positive constants, expres- 
sion (|3^|) reads 



1. (1-a) 1 

mr»= -lnw + 7 1 — 

8 a a 2 



<f 



(5) 



*|I 



(33) 



In obtaining the latter formula we have taken into account 
the equality 



ln(fc) e k dk = -7 . 



A similar expression can be obtained for the n-dimensional 
generalized Cauchy distribution [15] . 

Formula (|3"3")l immediately leads us to the desired ex- 
pression (fT4|) for W = ^mint- 



B Asymptotics of the probability density and 
generating function 

As in Appendix El let us consider a random vector r £ M. n 
whose statistical properties are specified by a given distri- 
bution density P(r). Now, however, some additional as- 
sumptions about it are adopted. Namely, the leading term 
in the asymptotic series of the probability density P(r) 
as r — > 00 is considered to be a spherically symmetrical 
power-law function. In other words, the equality 



P(r) 



j,n-\-a 



1 



j,n-\-a 



(34) 



where r = |r| and the parameters vo > and < a < 2 is 
presumed to hold. The inequality < a < 2 imposed on 
the exponent a assures us of the required convergence of 
the integral 

j P(r) dr = 1 

and meets the divergence of the second moment 
r 2 P(r) dr = 00 



being the characteristic feature of the systems at hand. 
The purpose of the present Appendix is to relate the given 
asymptotics of the probability density to corresponding 
asymptotics of the generating function G(k) as k — > 0. 
By virtue of ([22]) and ([34]) we have 



dr as k ->• (35) 



1 - G(k) = f (1 - e ikr ) P(r) dr 



2w 



sin 



r n+a 
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because due to the divergence of second moments the re- 
gion with |r| ~ l/|k| contributes mainly to integral f|35|) . 
A direct calculation of the last integral gives rise to 
the following asymptotic series of the generating function 

G(k) = l-w|k| a + o(|k| Q ) , (36) 

where the coefficient u is related to the parameter w as 

2 a . /7ra\ „f_ a\ „(n + a 



7T2 



sin , 

4f+i V 2 



r 1 



r 



(37) 



Dealing with a multidimensional space (n > 1) we can 
consider also the random variable r = |r| whose probabil- 
ity density P n (r) is related to the probability density P(r) 
of the random vector r as 



Pn(r) 



I P{r)da, 



(38) 



where the surface integral runs over the sphere §„ !r = 
{r : |r| = r}. In particular, for the probability density P(r) 
with the asymptotic behavior (|34|l the asymptotic series 
of probability density P n (f) is 



Pn(r) 



(39) 



where the coefficient w n = S n w and the cofactor S n is 
given by expression (|29p. Using (j3"Tf we can also write 



to 



2 qH 



/7ra\ „ / a\ 
5ln (-)^( 1+ 2) 



r 



"i 



(40) 



The obtained expressions can be interpreted in a re- 
verse manner. Namely, if the generating function G(k) 
of the random vector r exhibits the asymptotic behav- 
ior (j36|) . then the probability density functions P(r) and 
P n (r) of the random vector r itself and the derived ran- 
dom variable r = |r|, respectively, admit asymptotics (J34I) 
and (|39[) with the corresponding coefficients zu and w n 
related to the parameter uj by formulae (j3"7) and (|4"0"|) . 

Let us now measure the spatial scales of R ra in units of 
the geometric mean r g of the random vector r considered 
in Appendix \^ and assume the generating function G(r) 
to admit the asymptotic behavior (|36|) . Then, according 
to the obtained results, the random variable p = r/r g 
is characterized by the probability density V(p) with the 
asymptotics 



Hp) 



p 1 



where w p — zu n /r" and, thus, 

2 Q+1 /7ra\ / a\ r 
^ = — Sm {-) r { 1+ 2)- 



oo . 



n + a 



(41) 



, n 
n 2 



x exp < j(a — 1) 



#1 



*(?) 



(42) 



by virtue of ([3^| and (|40|) . Exactly this expression was 
used in writing expression (|17p . 
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